A stochastic approach for co-evolution process of virus and human immune system

Infectious diseases have long been a shaping force in human history, necessitating a comprehensive understanding of their dynamics. This study introduces a co-evolution model that integrates both epidemiological and evolutionary dynamics. Utilizing a system of differential equations, the model represents the interactions among susceptible, infected, and recovered populations for both ancestral and evolved viral strains. Methodologically rigorous, the model’s existence and uniqueness have been verified, and it accommodates both deterministic and stochastic cases. A myriad of graphical techniques have been employed to elucidate the model’s dynamics. Beyond its theoretical contributions, this model serves as a critical instrument for public health strategy, particularly predicting future outbreaks in scenarios where viral mutations compromise existing interventions.


Motivation
Our study advances the understanding of multiple infection dynamics, building on significant previous research in disease modeling.We resonate with and expand upon the work of researchers like those cited in 13,14 , and 15 , who explored transmission dynamics and sensitivity analysis using fractal-fractional differential operators.Our methodology also draws from 16,17 , who developed fractional models for diseases such as malaria and typhoid fever, underscoring the critical role of analyzing multiple pathogens in epidemiological studies and the importance of understanding their reproduction numbers.Furthermore, we recognize the impact of stochastic components in modeling diseases, as highlighted in 18 .The incorporation of time delays in stochastic epidemic models introduces complex system behaviors.The research in 19,20 , and 21 demonstrates how mathematical and statistical approaches can be employed to tackle questions of stochasticity and stability in epidemic models.Additionally, 22 's work on modeling COVID-19 with fractional order calculus emphasizes the effect of policy measures like isolation and vaccination on disease control, complementing our discussions on the efficacy of public health interventions.
Li and colleagues' contributions 23,24,25 shed light on the dynamic behavior of epidemic models, especially regarding bifurcations and chaotic phenomena, offering a valuable backdrop to our findings on oscillatory behaviors near endemic equilibrium.

Contribution
This work introduces a co-evolution model that seeks to capture this intricate balace between evolving viruses and the evolving human immune response.This model, built upon differential equations, represents the interaction dynamics of susceptible, infected, and recovered populations with both the original and evolved virus strains.By integrating evolutionary considerations into traditional epidemiological models, we aim to provide a more nuanced and comprehensive understanding of disease dynamics in the face of viral mutations and changing host immunities.

Novelty
Stochastic modeling plays a crucial role in epidemiology for several reasons.It allows researchers to account for the inherent randomness and variability in disease transmission among individuals and communities, reflecting the real-world unpredictability of outbreaks.This is vital for accurately simulating the spread of diseases and assessing potential outcomes under different scenarios.Stochastic models help in estimating the probabilities of different epidemiological events, enabling public health officials to make informed decisions regarding intervention strategies and resource allocation.Through stochastic modeling, epidemiologists can better understand the dynamics of infectious diseases, including the impact of factors like vaccination rates, population density, and social behaviors, thereby improving the effectiveness of disease control and prevention measures.
Beyond its theoretical significance, this model serves as a tool to guide public health interventions, and predict potential future outbreaks, especially in scenarios where rapid viral mutations might render existing treatments or vaccines less effective over time.In the subsequent sections, we will delve into the assumptions underlying the model, the mathematical formulations representing the co-evolution dynamics, and potential applications and implications of the findings derived from this model.

Model explanation
The model captures the interplay between human immunity and viral evolution.As the virus spreads, the human immune response evolves, leading to a changing transmission rate for the virus.The model can be used to analyze the impact of various interventions.

with
• S1 and S2 represent two subpopulations of the host species that are susceptible to infection by the virus strains Ȋ1 and Ȋ2 , respectively.• Ȓ represents the recovered individuals.
• β 1 and β 2 are the transmission rates of the two virus strains • r influences the rate at which the transmission rates ( β 1 and β 2 ) change in response to changes in the total number of infected individuals.• The parameter r significantly influences the transmission rates β 1 and β 2 for the two virus strains.Specifi- cally, the transmission rates are defined as indicating that the value of r directly affects how the total number of infected individuals I influences β 1 and β 2 .An increase in r amplifies the effect of changes in I on the transmission rates, due to the term α × r × I × 1 − I K .Thus, r acts as a modulator of the transmission rates' sensitivity to the infected population size, playing a crucial role in the disease's spread dynamics within the host population.
• The parameter α plays a crucial role in the model, acting as a factor within the feedback mechanism that influences the transmission rates β 1 and β 2 .Its primary function is to modulate the effect of the infected population size on the transmission rates, enabling the model to incorporate dynamic changes in transmission potential that occur as the prevalence of infection within the population changes.The model considers a moderate level of feedback, wherein increases in the infected population proportionally adjust the transmission rates, but not to an extreme degree.This adjustment is critical for accurately modeling the spread of infectious diseases, as it reflects the complex interactions between host behavior, population density, and pathogen transmissibility that can affect the rate at which an infection spreads through a community.• δ represents the death rate, µ is birth rate.
• ρ represents the rate at which recovered individuals lose immunity and move back into the susceptible cat- egory.
This models a loss of immunity over time.This work took motivation from recent work on stochastic modeling and infectious diseases models as discussed by [26][27][28][29][30][31][32] .The stochastic model is given by The model includes several σ parameters representing variability and stochastic effects: • σ 1 : Represents the variability in the susceptible population S1 , which arise from fluctuating contact rates or changes in population behavior that affect exposure to the first virus strain.• σ 2 : Captures the randomness in the second susceptible population S2 , which is due to similar factors as σ 1 , but with different underlying causes or magnitudes, given that S 2 represent a different risk group.• σ 3 : Reflects the random fluctuations in the number of individuals infected with the first virus strain Ȋ1 , due to variations in the disease's infectiousness, reporting rates, or response to treatment.• σ 4 : Pertains to the variability in the infection rate of the second virus strain Ȋ2 , which differ from σ 3 as the new strain has distinct characteristics, that is, higher transmissibility. ( (2) • σ 5 : Represents stochastic factors affecting the recovered population Ȓ , such as differential rates of loss of immunity or the impact of interventions that are not consistent across the entire population.
Each σ parameter is paired with a corresponding Brownian motion term, which mathematically represents the random 'noise' contributing to the fluctuations in each compartment over time.By including these stochastic terms, the model becomes a set of stochastic differential equations (SDEs), providing a more nuanced and realistic simulation of the epidemiological process, which can now capture both the average trends and the variability around those averages.The following assumptions underlie the model: • Every parameter within the system is a nonnegative, positive real number.
• The transmission rate of the virus can change based on the proportion of the infected population.
• Immunity to one strain does not confer immunity to the evolved strain.

Qualitative analysis
Definition 1 A system of differential equations is said to be Lipschitz continuous with respect to a variable x if there exists a constant L such that for any state variables x and y in a domain D , then , of the stochastic system.Then, for each equation in the system, the difference between the rates of change for the two solutions is bounded by a constant times the difference between the solutions.
Proof Taking equation for S1 , We can set, as an upper bound.
Repeating this process for the other equations, we can identify similar constants for each one.The largest of these constants then serves as the Lipschitz constant L for the whole system.Definition 2 A system is considered bounded if, for all solutions x(t) of the system and some positive constant M , the inequality ||x(t)|| ≤ M holds for all t.
Theorem 1 The stochastic system, given appropriate initial conditions, exhibits bounded behavior.

Proof Consider the first equation for S1 :
Given that populations cannot be negative, the loss terms ( β 1 S1 Ȋ1 , β 2 S1 Ȋ2 , δ S1 ) ensure that S1 does not grow unbounded.Therefore, the deterministic part of the system is bounded.

Theorem 2 Given boundedness and Lipschitz continuity, there exists a unique solution to the stochastic model for all time.
Proof For a system of stochastic differential equations (SDEs) of the form: where dW(t) represents the Wiener process, the existence and uniqueness of its solution is ensured if: • The coefficients f and g are bounded.
• The system is Lipschitz continuous.
We've shown that the system is Lipschitz continuous and bounded.Hence, we conclude that there exists a unique solution to the SDEs of our stochastic model for all time.

Equilibrium analysis
The variational matrix is obtained by linearizing the system of differential equations around the equilibrium.If we denote the equilibria as ( S 1e , S 2e , Ȋ 1e , Ȋ 2e , Re ) , the variational matrix V at this equilibrium is given by the Jacobian: E 1 (0, 0, 0, 0, 0) Given the variational matrix at the equilibrium E 1 (0, 0, 0, 0) , we have The eigenvalues of the variational matrix V at the equilibrium point E 1 (0, 0, 0, 0, 0 Therefore, the equilibrium E 1 (0, 0, 0, 0, 0) is semi-stable based on the eigenvalues of the vari- ational matrix.

Theorem 4
The stability point E 2 (1, 0, 0, 0, 0) exhibits local asymptotic stability provided that The remaining part of equilibrium analysis can be found in Appendix section 12.1.

Endemic equilibrium
The asymptotic solution relies heavily on the basic reproduction number R 0 of the disease, which is the expected number of cases directly generated by one case in a population where all individuals are susceptible to infection 33 .The presence of an endemic equilibrium depends on various factors, including the basic reproduction number R 0 .The Jacobian matrix J of the system at the DFE (disease free equilibrium) is given by: The next-generation matrix, K , is given by the product of two matrices, F and V −1 , For our system, the matrix F is . The inverse is The next-generation matrix K is Thus, R 0 is the maximum of the two diagonal entries.
For R 01 for Ȋ1 is β 2 S1 γ +σ +δ , For R 02 for Ȋ2 is β 1 S2 γ +σ +δ .To determine the endemic equilibrium, we will evaluate the following basic reproduction numbers, • R I 1 represents the mean count of recent infection cases attributed to I 1 .
• R I 2 symbolizes the average count of new infection cases ascribed to I 2 .
The system's critical parameter can be expressed as, Beyond the equilibria highlighted earlier, the model's endemic equilibrium is realized when, Theorem 5 There exists a unique endemic equilibrium whenever R 0 > 1.
Proof can be found in Appendix section 12.2.

Stochastic analysis
Suppose a probability domain represented as (�, G, Q) containing a Wiener process (or Brownian motion) rep- resented as W = W π , G W π , π > 0 .The associated filtration is given by (G π , π > 0) .Let as the governing stochastic differential equation.The function v(π, Z(π)) maps from [0, ∞) × R e to R e .y(π, Z(π)) is considered to be an p × q matrix.In this scenario, both y and v satisfy Lipschitz conditions.Let's introduce K as the differential operator for the system described in equation 5 , Applying operator K on a function , where Our aim is to determine κ 0 and κ > 0 such that, t = 0 a.s., we obtain, For the stochastic system with S1 , S2 , Ȋ1 , Ȋ2 , R , we define as, where X i is the corresponding compartment.It remains to prove that adheres to the a.s., invariance principle.Proof can be found in Appendix section 12.3.
Theorem 7 adheres to the almost sure invariance principle for model (1).
Proof can be found in Appendix section 12.4.
Proof can be found in Appendix section 12.5.

Analysis of the population dynamics
The following values in Table 1 has been used for numerical simulation.These values are based on theoretical studies and empirical findings [35][36][37][38]

Time series
The Euler-Maruyama method for a stochastic differential equation (SDE) of the form is given by the following recursive formula, where a(X t , t) is the drift coefficient, b(X t , t) is the diffusion coefficient, and W n is the Wiener process incre- ment approximated by √ �t • N(0, 1) .The Euler-Maruyama approximation of the given stochastic model is, • The susceptible populations, represented by S1 and S2 , manifest a declining trend in Fig. 3.This reduction is indicative of the individuals transitioning to the infected categories, depleting the pool of individuals that can potentially be infected.• The infected categories, represented by Ȋ1 and Ȋ2 , show a typical infectious disease trajectory, an initial rise as the disease spreads through the susceptible population, a peak representing the maximum number of concurrent infections, and a subsequent decline as individuals recover or die.• The Ȓ category, which denotes recovered individuals, exhibits an increasing trend, reflecting the accumula- tion of individuals who have gained immunity post-infection.
In Fig. 2, the deterministic lines depict an expected trajectory of populations, revealing trends of decline in susceptible classes, a peak in infections, and a rise in recoveries.In contrast, the stochastic lines mirror these general trends but with evident fluctuations, representing real-world variability due to unpredictable factors.Specifically, while susceptibles for both ancestral and evolved strains decrease over time, the evolved strain depletes the susceptible pool more rapidly.The number of infections shows varied dynamics, with the ancestral strain presenting a rapid rise and decline and the evolved strain maintaining a more consistent infection rate.The recovered population steadily increases over time.The disparities between deterministic and stochastic representations underscore the importance of factoring in real-world randomness alongside average predictions in infectious disease modeling.
Table 1.Initial and parametric values.www.nature.com/scientificreports/Ȋ1 is less virulent or less successful at sustaining its spread compared to Ȋ2 given that Ȋ1 dies out, while Ȋ2 reaches a steady state.S1 and S2 both experience significant drops, but S1 has a slight recovery, showing resistance in the population against Ȋ1 .The Ȓ compartment's steady rise and plateau suggest that recovery is the primary outcome for infected individuals, which is a positive sign in terms of public health.2, drawing on extensive datasets described in contemporary literature 11,35,[39][40][41][42][43] .
We have employed the nonlinear least square method for the estimation process.The Ordinary Least Square (OLS) solution was implemented to minimize the error terms as delineated in Eq. ( 7), and the related relative error was employed to assess the goodness of fit. Figure 3 illustrates the infected population as predicted by our proposed system with the real represented by stars.Here the notion I i is the reported cumulative infected cases and Îi is the cumulative infected cases obtained from simulating the model.

Phase plane plot
The trajectories in Fig. 4 reveal the cyclical nature of the interactions between susceptibility and infection, illuminating the ebb and flow of the epidemic.As susceptibles are depleted, the number of new infections decreases, leading to a decline in the infected populations.While Ȋ1 have a cyclical pattern with S1 , indicating of outbreak and control, Ȋ2 shows a more direct path to depletion of susceptibles, suggesting a more virulent or transmissible strain.This underlines the importance of strain-specific public health interventions.

Time series with moving average
While the raw data for Ȋ1 and Ȋ2 captures the actual number of infections, the moving average smoothens out short-term fluctuations, thereby highlighting broader trends.Both the actual infection data and the moving average exhibit a prominent peak, representing the height of the epidemic before a decline sets in.This underscores the transient nature of outbreaks and the eventual return to equilibrium as shown by Fig. 5.

Quiver plot
The quiver plot in Fig. 6 suggests complex interactions between the Ȋ1 and Ȋ2 strains.The arrows elucidate the direction and magnitude of change, providing a snapshot of the evolution in infection rates.The looping trajectory show oscillations in the infected populations, suggesting periodic outbreaks or potential for repeated waves of infections.This cyclic pattern arises from various factors, such as loss of immunity, seasonal variations, or reintroduction of infections.

Population distribution
This bar graph gives an overview of the distribution of different populations at a specific time, t = 100 in Fig. 7.By this time, a significant portion of the initial susceptible population have transitioned to the recovered category.The number of infected individuals ( Ȋ1 and Ȋ2 ) is notably less than the susceptible and recovered groups, indicating that the epidemic is waning.

Conclusion
The introduced epidemiological model represents a significant enhancement in the mathematical study of infectious diseases, especially those with multiple circulating strains such as SARS-CoV-2.Authenticated by real data analysis, the model's equations provide a precise depiction of the disease dynamics, offering authenticity and bolstering its credibility within the scientific community.The analysis of Ȋ1 and Ȋ2 through this model showcases a comprehensive mathematical characterization of the spread and potential control of multiple viral strains.The incorporation of empirical data not only underscores the model's accuracy but also fortifies its predictive capacity regarding disease progression and potential  extinction scenarios.The model's robustness is affirmed through rigorous exploration of disease extinction conditions, laying a foundation for understanding the critical thresholds that govern the eradication of the disease.
The research by 10 models COVID-19's transmission with a focus on variant dynamics and vaccine impact, while 11 assesses the influence of variants in France using optimization for parameter estimation.The study in 12 introduces a fractional model for two COVID-19 strains, and 1 explores the impact of multiple strains on pandemic trajectories and vaccine efficacy.These studies contribute valuable insights into pandemic modeling; however, our work extends these efforts by offering a more intricate examination of disease dynamics, specially; • Integrates a wider range of epidemiological and evolutionary dynamics.
• Provides a deeper analytical approach including a comprehensive equilibrium analysis, and disease extinction conditions.• Offers a more detailed exploration of both deterministic and stochastic scenarios.
• Utilizes advanced graphical techniques for a clearer understanding of disease progression.
• Addresses practical applications for real-world outbreaks, especially in the context of evolving viral strains impacting public health measures.

Remarks and future recommendations
Our in-depth examination of the co-evolution model of host and human immune response highlights the intricate dynamics between infectious strains Ȋ1 and Ȋ2 .The stark differences in the trajectories of these strains underscore the complex challenges posed in managing multi-strain infectious diseases.While our model offers a comprehensive understanding, the continually evolving nature of infectious diseases calls for adaptive strategies and persistent refinement in modeling approaches.
For future studies, it is recommended to incorporate factors like varying mutation rates, potential crossimmunity effects, and the impact of external interventions such as vaccination campaigns.Additionally, integrating real-world data can augment the model's predictive accuracy.As global communities grapple with the multifaceted challenges posed by infectious diseases, models such as ours serve as foundational tools, and their continuous refinement remains pivotal for informed public health decision-making.

Figure 1 .
Figure 1.Virus and immune system's co-evolution process.

Figure 5 .
Figure 5.Time series with moving average.